This notebook uses changes in the Normalized Difference Vegetation Index (NDVI) to identify vegetation change. The algorithm identifies a "baseline" and "analysis" time period and then compares the spectral index in each of those time periods. Significant changes in NDVI (vegetation greenness) are coincident with land change, as long as the comparisons are done between similar time periods (seasons or years). Users of this algorithm should not accept the accuracy of the results but should conduct ground validation testing to assess accuracy. It is expected that this algorithm can be used to identify clusters of pixels that have experienced change and allow targeted investigation of those areas by local or regional governments. In some cases the impacts may be negative (deforestation, mining, burning, drought) or positive (regrowth, improved soil moisture).
It should also be noted that the selection of the baseline and analysis time period is critical. First, the two time periods should be similar (season, year) so that the vegetation state can be compared in similar weather conditions. Second, the time periods should be sufficiently clear (non-cloudy) data. If the baseline or analysis mosaic (composite of images) is contaminated with clouds, it will impact the results.
import datacube
dc = datacube.Datacube(app='Vegetation_Change')
import sys, os
os.environ['USE_PYGEOS'] = '0'
from datacube.utils import masking
from dea_tools.plotting import rgb, display_map
from odc.algo import to_f32
### EASI tools
sys.path.append(os.path.expanduser('../scripts'))
from ceos_utils.data_cube_utilities.clean_mask import landsat_clean_mask_invalid, landsat_qa_clean_mask
from easi_tools import EasiDefaults
from easi_tools import notebook_utils
easi = EasiDefaults() # Get the default parameters for this system
Successfully found configuration for deployment "eail"
cluster, client = notebook_utils.initialize_dask(use_gateway=False)
display(cluster if cluster else client)
print(notebook_utils.localcluster_dashboard(client, server=easi.hub))
64a0b9b9
| Dashboard: http://127.0.0.1:8787/status | Workers: 4 |
| Total threads: 8 | Total memory: 29.00 GiB |
| Status: running | Using processes: True |
Scheduler-1d3a4938-5207-40b6-aa9c-3b0bbb6c36c3
| Comm: tcp://127.0.0.1:43369 | Workers: 4 |
| Dashboard: http://127.0.0.1:8787/status | Total threads: 8 |
| Started: Just now | Total memory: 29.00 GiB |
| Comm: tcp://127.0.0.1:42903 | Total threads: 2 |
| Dashboard: http://127.0.0.1:39163/status | Memory: 7.25 GiB |
| Nanny: tcp://127.0.0.1:40295 | |
| Local directory: /tmp/dask-worker-space/worker-twl206p5 | |
| Comm: tcp://127.0.0.1:37309 | Total threads: 2 |
| Dashboard: http://127.0.0.1:42975/status | Memory: 7.25 GiB |
| Nanny: tcp://127.0.0.1:42435 | |
| Local directory: /tmp/dask-worker-space/worker-7cansc5h | |
| Comm: tcp://127.0.0.1:38089 | Total threads: 2 |
| Dashboard: http://127.0.0.1:38395/status | Memory: 7.25 GiB |
| Nanny: tcp://127.0.0.1:44559 | |
| Local directory: /tmp/dask-worker-space/worker-pz6usymu | |
| Comm: tcp://127.0.0.1:37713 | Total threads: 2 |
| Dashboard: http://127.0.0.1:37015/status | Memory: 7.25 GiB |
| Nanny: tcp://127.0.0.1:46389 | |
| Local directory: /tmp/dask-worker-space/worker-damtr0hx | |
https://hub.eail.easi-eo.solutions/user/jhodge/proxy/8787/status
from datacube.utils.aws import configure_s3_access
configure_s3_access(aws_unsigned=False, requester_pays=True, client=client)
<botocore.credentials.Credentials at 0x7f048557d780>
# Select a Product
product = "landsat8_c2l2_sr"
# Select an analysis region (Latitude-Longitude)
# Mining Region near Obuasi, Ghana
# Baseline = 12/23/13, Analysis = 1/6/19
# NOTE: These are clear scenes at similar times of the year.
# Use the Cloud Statistics notebook to evaluate cloud cover.
latitude = (-23.5,-23.6)
longitude = (-59.45,-59.35)
## The code below renders a map that can be used to orient yourself with the region.
display_map(longitude,latitude)
/env/lib/python3.10/site-packages/dea_tools/plotting.py:313: FutureWarning: This function is deprecated. See: https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1
all_longitude, all_latitude = transform(Proj(crs), Proj('EPSG:4326'), all_x,
from datetime import datetime
# Select the start and end periods for your analysis products (Year,Month,Day).
# These time windows will be used to make a median mosaic, so typically pick a year length (or more)
# or select a small window surrounding a clear single date (use Cloud Statistics notebook).
# Also, be sure to evaluate the RGB mosaics (below) to affirm they are not full of clouds.
# Select the baseline time period (start and end)
baseline_time_period = (datetime(2013,11,12), datetime(2013,11,13))
# Select the analysis time period (start and end)
analysis_time_period = (datetime(2022,12,23), datetime(2022,12,24))
output_crs = 'epsg:6933'
resolution = (-30,30)
scale = 0.0000275
offset = 0.2
baseline_ds = dc.load(latitude=latitude,longitude=longitude,product=product,
measurements = ['red', 'green', 'blue', 'nir', 'swir1', 'swir2', 'pixel_qa'],
time=baseline_time_period,output_crs=output_crs,resolution=resolution,group_by='solar_day',dask_chunks={'time':1})
analysis_ds = dc.load(latitude=latitude,longitude=longitude,product=product,
measurements = ['red', 'green', 'blue', 'nir', 'swir1', 'swir2', 'pixel_qa'],
time=analysis_time_period,output_crs=output_crs,resolution=resolution,group_by='solar_day',dask_chunks={'time':1})
cloud_mask_baseline = masking.make_mask(baseline_ds['pixel_qa'], clear='clear', cloud='not_high_confidence', cloud_shadow='not_high_confidence')
baseline_ds = baseline_ds.where(cloud_mask_baseline)
cloud_mask_analysis = masking.make_mask(analysis_ds['pixel_qa'], cloud='not_high_confidence', cloud_shadow='not_high_confidence')
analysis_ds = analysis_ds.where(cloud_mask_analysis)
from ceos_utils.data_cube_utilities.dc_mosaic import create_median_mosaic
baseline_composite = to_f32(create_median_mosaic(baseline_ds, cloud_mask_baseline), scale=scale, offset=offset)
analysis_composite = to_f32(create_median_mosaic(analysis_ds, cloud_mask_analysis), scale=scale, offset=offset)
def NDVI(dataset):
return (dataset.nir - dataset.red)/(dataset.nir + dataset.red)
parameter_baseline_composite = NDVI(baseline_composite)
parameter_analysis_composite = NDVI(analysis_composite)
parameter_anomaly = parameter_analysis_composite - parameter_baseline_composite
import matplotlib.pyplot as plt
from ceos_utils.data_cube_utilities.dc_rgb import rgb
from matplotlib.cm import RdYlGn
RdYlGn.set_bad('black',1.)
# Define the significant anomaly range for Plot #4
# The typical loss range is <0. Choose < -0.1 or < -0.2.
# The typical gain range is >0. Choose > 0.1 or > 0.2.
loss_range = parameter_anomaly < -0.05
gain_range = parameter_anomaly > 0.1
import xarray as xr
import numpy as np
fig, ax = plt.subplots(2, 2, figsize=(12,12))
for sub_ax in ax.flatten():
sub_ax.set_facecolor('black')
baseline_rgb = baseline_composite[['red', 'green', 'blue']].to_array().compute()
analysis_rgb = analysis_composite[['red', 'green', 'blue']].to_array().compute()
# Use the middle values of the data (2% to 98%) to brighten the image
lw_qtl, up_qtl = 0.02, 0.98
rgb_vmin = min(baseline_rgb.quantile(lw_qtl).values,analysis_rgb.quantile(lw_qtl).values)
rgb_vmax = max(baseline_rgb.quantile(up_qtl).values,analysis_rgb.quantile(up_qtl).values)
# Plot the resulting 4 products ... Baseline RGB, Analysis RGB, Total Anomaly, Anomaly Threshold
# NOTE: Clouds in either the baseline or analysis images will be removed from the anomaly product
## Plot #1 = Baseline RGB (upper left)
axes_image = baseline_rgb.plot.imshow(ax=ax[0,0], vmin=rgb_vmin, vmax=rgb_vmax)
## Plot #2 = Analysis RGB (upper right)
analysis_rgb.plot.imshow(ax=ax[0,1], vmin=rgb_vmin, vmax=rgb_vmax)
## Plot #3 = Total Anomaly (lower left)
parameter_anomaly.plot(ax=ax[1,0], vmin=-0.2, vmax=0.2, cmap = RdYlGn, add_colorbar=False)
## Plot #4 = Anomaly Threshold (lower right)
# Analysis composite grayscale background
plt4_bkg_band = 'swir1' # The band to use as the background image.
plt4_rgb = np.repeat(analysis_composite[plt4_bkg_band].where(cloud_mask_baseline.sum('time').astype(np.bool))\
.values[:,:,np.newaxis],3,axis=2)
# Selected a range of SWIR1 values (0.001 to 0.600) to lighten image background
# Users may also try values of 0.02 and 0.98
min_bkg = np.nanquantile(analysis_composite[plt4_bkg_band].values, 0.001)
max_bkg = np.nanquantile(analysis_composite[plt4_bkg_band].values, 0.95)
plt4_rgb = np.interp(plt4_rgb, (min_bkg, max_bkg), [0,1])
# Significant anomaly color overlays
color_green = np.array([0,1,0]) # green
color_red = np.array([1,0,0]) # red
plt4_rgb[loss_range] = color_red
plt4_rgb[gain_range] = color_green
# Plot
plt4_coords = dict(analysis_composite.coords)
rgb_coord_arr = np.array(['red', 'green', 'blue'])
rgb_coord_da = xr.DataArray(rgb_coord_arr,name='rgb',dims=['rgb'],coords={'rgb': rgb_coord_arr})
plt4_coords.update({'rgb': rgb_coord_da})
plt4_rgb_da = xr.DataArray(plt4_rgb, coords=plt4_coords,dims=list(analysis_composite.dims) + ['rgb'])
plt4_rgb_da.plot.imshow(ax=ax[1,1])
# Titles for all plots
ax[0,0].set_title('Baseline Composite'), ax[0,0].xaxis.set_visible(False), ax[0,0].yaxis.set_visible(False)
ax[0,1].set_title('Analysis Composite'), ax[0,1].xaxis.set_visible(False), ax[0,1].yaxis.set_visible(False)
ax[1,0].set_title('Vegetation Anomalies: Red=Loss, Green=Gain'), ax[1,0].xaxis.set_visible(False), ax[1,0].yaxis.set_visible(False)
ax[1,1].set_title('Locations of Significant Anomalies: Red=Loss, Green=Gain'), ax[1,1].xaxis.set_visible(False), ax[1,1].yaxis.set_visible(False)
plt.tight_layout()
plt.show()
/tmp/ipykernel_510/2898797902.py:33: DeprecationWarning: `np.bool` is a deprecated alias for the builtin `bool`. To silence this warning, use `bool` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.bool_` here.
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
plt4_rgb = np.repeat(analysis_composite[plt4_bkg_band].where(cloud_mask_baseline.sum('time').astype(np.bool))\